City-level information pulled from the following sources:
LA Almanac (race, age, and income) http://www.laalmanac.com/population/po38.php http://www.laalmanac.com/employment/em12.php
LA County vaccine percent by city http://publichealth.lacounty.gov/media/coronavirus/vaccine/vaccine-dashboard.htm
LA County age adjusted incidence rates http://dashboard.publichealth.lacounty.gov/covid19_surveillance_dashboard/
Data are incidence rates (proportions), vaccination rates (proportions), and characteristics by city.
All hypothesis are on the city-level and care must be taken when interpreting the results as they describe city-level variation and not individual level variation or co-variation.
Investigation of the December 2020-January 2021 peak for COVID 19 captured by a single day for incidence rates: 2021-01-07:
Investigation of the August 2021 peak for COVID 19 captured by a single day for incidence rates: 2021-08-18:
Investigation of August 2021 vaccination rates for COVID 19 captured by a single day for vaccination rates: 2021-08-18:
d <- readRDS("data/analysis_data.rds") %>%
arrange(Date.x) %>%
filter(Date.x == "2021-08-18")
city.populations <- data.frame(d$city, d$population.x)
names(city.populations) <- c("City", "Population_Size")
city.populations <- city.populations[order(city.populations$Population_Size, decreasing=T), ]
kable(city.populations, format = "simple", align = 'c', caption = "Cities and Population Size")
| City | Population_Size | |
|---|---|---|
| 121 | Los Angeles | 4043337 |
| 95 | Santa Clarita | 220424 |
| 44 | Glendale | 206493 |
| 62 | Lancaster | 161570 |
| 78 | Palmdale | 158968 |
| 82 | Pomona | 157869 |
| 108 | Torrance | 149268 |
| 35 | East Los Angeles | 125269 |
| 39 | El Monte | 117269 |
| 33 | Downey | 114263 |
| 52 | Inglewood | 113582 |
| 114 | West Covina | 108235 |
| 77 | Norwalk | 107622 |
| 18 | Burbank | 107180 |
| 25 | Compton | 99904 |
| 101 | South Gate | 98155 |
| 20 | Carson | 93846 |
| 97 | Santa Monica | 92446 |
| 48 | Hawthorne | 91301 |
| 119 | Whittier | 91216 |
| 4 | Alhambra | 86724 |
| 61 | Lakewood | 80362 |
| 15 | Bellflower | 77735 |
| 12 | Baldwin Park | 76769 |
| 68 | Lynwood | 72047 |
| 85 | Redondo Beach | 68697 |
| 11 | Azusa | 65963 |
| 26 | Covina | 65851 |
| 6 | Arcadia | 65735 |
| 42 | Florence-Graham | 64705 |
| 74 | Montebello | 64375 |
| 81 | Pico Rivera | 64284 |
| 75 | Monterey Park | 62262 |
| 43 | Gardena | 61310 |
| 51 | Huntington Park | 59484 |
| 104 | South Whittier | 59222 |
| 32 | Diamond Bar | 57515 |
| 80 | Paramount | 56023 |
| 46 | Hacienda Heights | 55926 |
| 88 | Rosemead | 55350 |
| 45 | Glendora | 52764 |
| 89 | Rowland Heights | 51022 |
| 22 | Cerritos | 50067 |
| 56 | La Mirada | 49599 |
| 5 | Altadena | 43620 |
| 14 | Bell Gardens | 43071 |
| 84 | Rancho Palos Verdes | 42747 |
| 73 | Monrovia | 42681 |
| 8 | Westmont | 42442 |
| 92 | San Gabriel | 40954 |
| 57 | La Puente | 40697 |
| 29 | Culver City | 39865 |
| 115 | West Hollywood | 36951 |
| 23 | Claremont | 36484 |
| 107 | Temple City | 36455 |
| 13 | Bell | 36332 |
| 70 | Manhattan Beach | 35999 |
| 58 | La Verne | 35322 |
| 120 | Willowbrook | 34913 |
| 16 | Beverly Hills | 34520 |
| 90 | San Dimas | 34516 |
| 63 | Lawndale | 33614 |
| 111 | Walnut | 30532 |
| 72 | Maywood | 28049 |
| 21 | Castaic | 27191 |
| 34 | Duarte | 26444 |
| 102 | South Pasadena | 26053 |
| 91 | San Fernando | 24612 |
| 28 | Cudahy | 24347 |
| 19 | Calabasas | 24323 |
| 76 | East San Gabriel | 24036 |
| 110 | Valinda | 23371 |
| 100 | South El Monte | 22680 |
| 64 | Lennox | 22542 |
| 113 | West Carson | 22086 |
| 105 | Stevenson Ranch | 20966 |
| 2 | Agoura Hills | 20883 |
| 67 | Lomita | 20729 |
| 54 | La Crescenta-Montrose | 19801 |
| 49 | Hermosa Beach | 19670 |
| 96 | Santa Fe Springs | 18364 |
| 7 | Artesia | 16795 |
| 40 | El Segundo | 16786 |
| 112 | Walnut Park | 16143 |
| 37 | East Rancho Dominguez | 15308 |
| 47 | Hawaiian Gardens | 14676 |
| 79 | Palos Verdes Estates | 13522 |
| 93 | San Marino | 13277 |
| 27 | Charter Oak | 13144 |
| 24 | Commerce | 13069 |
| 60 | Lake Los Angeles | 12994 |
| 69 | Malibu | 12961 |
| 83 | Quartz Hill | 12906 |
| 99 | Signal Hill | 11797 |
| 98 | Sierra Madre | 10989 |
| 116 | West Puente Valley | 9835 |
| 71 | Marina del Rey | 9411 |
| 103 | South San Gabriel | 8848 |
| 118 | Westlake Village | 8360 |
| 87 | Rolling Hills Estates | 8113 |
| 1 | Acton | 7971 |
| 59 | Ladera Heights | 7071 |
| 10 | Avocado Heights | 6775 |
| 36 | East Pasadena | 6403 |
| 106 | Sun Village | 6036 |
| 55 | La Habra Heights | 5455 |
| 38 | East Whittier | 5306 |
| 30 | Del Aire | 4393 |
| 3 | Agua Dulce | 4158 |
| 66 | Littlerock | 4021 |
| 9 | Avalon | 3869 |
| 109 | Val Verde | 3309 |
| 31 | Desert View Highlands | 2493 |
| 94 | San Pasqual | 2035 |
| 86 | Rolling Hills | 1940 |
| 50 | Hidden Hills | 1890 |
| 65 | Leona Valley | 1751 |
| 41 | Elizabeth Lake | 1661 |
| 53 | Irwindale | 1459 |
| 117 | West Rancho Dominguez | 1359 |
| 17 | Bradbury | 1069 |
City <- d$city
Jan.IR <- d$adj_case_14day_rate.y
Jan.IR.cat <- relevel(as.factor(ifelse(Jan.IR > quantile(Jan.IR, probs=.75), "High",
ifelse(Jan.IR < quantile(Jan.IR, probs=.25), "Low",
"Med"))), ref="Med")
hist(Jan.IR, breaks=seq(from=0, to=round(max(Jan.IR),-3), by=200), col="red", main="January Incidence Rates for Cities", xlab="Incidence Rates")
abline(v=quantile(Jan.IR, probs=c(.25, .5, .75)), lwd=2, lty=2)
Aug.IR <- d$adj_case_14day_rate.x
hist(Aug.IR, breaks=seq(from=0, to=round(max(Aug.IR),-2), by=100), col="red", main="August Incidence Rates for Cities", xlab="Incidence Rates")
Aug.Vax.percent <- d$dose1_all_c_prcent
Aug.Vax.cat <- relevel(as.factor(ifelse(Aug.Vax.percent > quantile(Aug.Vax.percent, probs=.75), "High",
ifelse(Aug.Vax.percent < quantile(Aug.Vax.percent, probs=.25), "Low",
"Med"))), ref="Med")
hist(Aug.Vax.percent, breaks=seq(from=0, to=100, by=5), col="blue", main="Percent Vaccination Rates for Cities", xlab="Percent Vaccinated")
abline(v=quantile(Aug.Vax.percent, probs=c(.25, .5, .75)), lwd=2, lty=2)
CityHispanic <- d$Hispanic
CityHispanic.cat <- relevel(as.factor(ifelse(d$Hispanic > quantile(d$Hispanic, probs=.75), "High",
ifelse(d$Hispanic < quantile(d$Hispanic, probs=.25), "Low",
"Med"))), ref="Med")
hist(CityHispanic, breaks=seq(from=0, to=100, by=5), main="Proportion of Hispanic Individuals", xlab="Proportion")
abline(v=quantile(CityHispanic, probs=c(.25, .5, .75)), lwd=2, lty=2)
CityBlack <- d$Black
CityBlack.cat <- relevel(as.factor(ifelse(d$Black > 10, "High","Low")), "Low")
hist(CityBlack, breaks=seq(from=0, to=100, by=5), main="Proportion of Black Individuals", xlab="Proportion")
abline(v=10, lwd=2, lty=2)
text(10, 60, "Cutoff at 10%", pos=4)
CityAsian <- d$Asian
CityAsian.cat <- relevel(as.factor(ifelse(d$Asian > quantile(d$Asian, probs=.75), "High",
ifelse(d$Asian < quantile(d$Asian, probs=.5), "Low",
"Med"))), ref="Low")
hist(CityAsian, breaks=seq(from=0, to=100, by=5), main="Proportion of Asian Individuals", xlab="Proportion")
abline(v=quantile(CityAsian, probs=c(.25, .5, .75)), lwd=2, lty=2)
Age data is captured by the number/proportion of individuals within each city in a age category. Age categories include: “Age.Under.15”,“Age.15.17”,“Age.18.24”,“Age.25.34”,“Age.35.54”,“Age.55.64”,“Age.65.”
Three groups are created for the analysis:
Young: combined age categories for “Age.Under.15”,“Age.15.17”,“Age.18.24”
Middle: age category for “Age.35.54”
Old: combined age categories for “Age.55.64”,“Age.65.”
Age <- d[,c("Age.Under.15","Age.15.17","Age.18.24","Age.25.34","Age.35.54","Age.55.64","Age.65.")]
Age <- t(apply(Age, 1, FUN=function(v) { v/sum(v) }))
CityAgeYoung <- apply(Age[,c("Age.Under.15","Age.15.17","Age.18.24")], 1, sum)
CityAgeOld <- apply(Age[,c("Age.55.64","Age.65.")], 1, sum)
CityAgeMid <- Age[,"Age.35.54"]
CityAge <- cbind(CityAgeYoung, CityAgeOld)
CityAgeOld.cat <- relevel(as.factor(ifelse(CityAgeOld > quantile(CityAgeOld, probs=.75), "High","notHigh")), "notHigh")
CityAgeYoung.cat <- relevel(as.factor(ifelse(CityAgeYoung > quantile(CityAgeYoung, probs=.75), "High","notHigh")), "notHigh")
ftable(CityAgeOld.cat, CityAgeYoung.cat)
## CityAgeYoung.cat notHigh High
## CityAgeOld.cat
## notHigh 61 30
## High 30 0
HouseholdIncome <- d$Households
HouseholdIncome.cat <- relevel(as.factor(ifelse(HouseholdIncome > quantile(HouseholdIncome, probs=.75), "High",
ifelse(HouseholdIncome < quantile(HouseholdIncome, probs=.25), "Low",
"Med"))), ref="Med")
hist(HouseholdIncome, breaks=seq(from=0, to=round(max(HouseholdIncome),-2), by=10000), main="Household Income Within Each City", xlab="Income", col="green")
abline(v=quantile(HouseholdIncome, probs=c(.25, .5, .75)), lwd=2, lty=2)
CityAge.m <- cbind(CityAgeOld.cat, CityAgeYoung.cat)
CityRaceEthnicty.m <- cbind(CityHispanic.cat, CityBlack.cat, CityAsian.cat)
d.analysis <- data.frame(Aug.IR, Jan.IR.cat, Aug.Vax.cat, Aug.Vax.percent, CityRaceEthnicty.m, CityAge.m, HouseholdIncome.cat)
summarytools::view(dfSummary(d.analysis, style = 'grid',
max.distinct.values = 10, plain.ascii = FALSE, valid.col = FALSE, headings = FALSE), method = "render")
| No | Variable | Stats / Values | Freqs (% of Valid) | Graph | Missing | ||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | Aug.IR [numeric] | Mean (sd) : 436.5 (178.6) min < med < max: 39 < 447 < 1086 IQR (CV) : 186 (0.4) | 110 distinct values | 0 (0.0%) | |||||||||||||||||||
| 2 | Jan.IR.cat [factor] | 1. Med 2. High 3. Low |
|
0 (0.0%) | |||||||||||||||||||
| 3 | Aug.Vax.cat [factor] | 1. Med 2. High 3. Low |
|
0 (0.0%) | |||||||||||||||||||
| 4 | Aug.Vax.percent [numeric] | Mean (sd) : 63.7 (11.2) min < med < max: 19.3 < 64.9 < 83.8 IQR (CV) : 14.6 (0.2) | 104 distinct values | 0 (0.0%) | |||||||||||||||||||
| 5 | CityHispanic.cat [integer] | Mean (sd) : 1.7 (0.8) min < med < max: 1 < 1 < 3 IQR (CV) : 1 (0.5) |
|
0 (0.0%) | |||||||||||||||||||
| 6 | CityBlack.cat [integer] | Min : 1 Mean : 1.1 Max : 2 |
|
0 (0.0%) | |||||||||||||||||||
| 7 | CityAsian.cat [integer] | Mean (sd) : 1.8 (0.8) min < med < max: 1 < 2 < 3 IQR (CV) : 2 (0.5) |
|
0 (0.0%) | |||||||||||||||||||
| 8 | CityAgeOld.cat [integer] | Min : 1 Mean : 1.2 Max : 2 |
|
0 (0.0%) | |||||||||||||||||||
| 9 | CityAgeYoung.cat [integer] | Min : 1 Mean : 1.2 Max : 2 |
|
0 (0.0%) | |||||||||||||||||||
| 10 | HouseholdIncome.cat [factor] | 1. Med 2. High 3. Low |
|
0 (0.0%) |
Generated by summarytools 0.9.9 (R version 4.1.0)
2021-11-04
# need to make sure low is lowest when setting as.numeric (to make cor more clear)
d.num <- data.frame(Aug.IR, Jan.IR.cat, Aug.Vax.cat,
CityHispanic.cat, as.numeric(CityBlack.cat), as.numeric(CityAsian.cat),
as.numeric(CityAgeOld.cat), as.numeric(CityAgeYoung.cat),
HouseholdIncome.cat) %>%
mutate(across(c(Jan.IR.cat, Aug.Vax.cat, CityHispanic.cat, HouseholdIncome.cat),
~ as.numeric(relevel(.x, ref = "Low"))))
cormat <- cor(d.num, use="complete.obs", method="pearson")
corrplot(cormat, type="upper", order="original",
col=brewer.pal(n=8, name="RdYlBu"),
title = "",
addCoef.col = "black",
tl.cex=.5, number.cex=.5)
CityHispanic.cont <- as.numeric(relevel(CityHispanic.cat, "Low"))
coef <- summary(lm(Jan.IR~CityHispanic.cont))$coef
ggplot(d.analysis, aes(x=CityHispanic.cont, y=Jan.IR)) +
geom_point(color="blue", size=2) +
geom_smooth(method=lm, color="black") +
geom_text(x=2.5, y=3000, label=paste(expression(beta), "=", round(coef[2,1],2))) +
geom_text(x=2.5, y=2500, label=paste("p-value", "=", round(coef[2,4],25))) +
scale_x_continuous(name ="Proportion of Hispanic Individuals", breaks=c(1,2,3), labels=c("Low", "Med", "High"))+
ylab("Jan Incidence Rates")
## `geom_smooth()` using formula 'y ~ x'
coef <- summary(lm(Jan.IR~as.numeric(relevel(HouseholdIncome.cat, "Low"))))$coef
ggplot(d.analysis, aes(x=as.numeric(relevel(HouseholdIncome.cat, "Low")), y=Jan.IR)) +
geom_point(color="blue", size=2) +
geom_smooth(method=lm, color="black") +
geom_text(x=2.5, y=3000, label=paste(expression(beta), "=", round(coef[2,1],2))) +
geom_text(x=2.5, y=2500, label=paste("p-value", "=", round(coef[2,4],23))) +
scale_x_continuous(name ="Household Income", breaks=c(1,2,3), labels=c("Low", "Med", "High"))+
ylab("January Incidence Rates")
## `geom_smooth()` using formula 'y ~ x'
coef <- summary(lm(Aug.IR~Jan.IR))$coef
ggplot(d.analysis, aes(x=Jan.IR, y=Aug.IR)) +
geom_point(color="blue", size=2) +
geom_smooth(method=lm, color="black") +
geom_text(x=3000, y=900, label=paste(expression(beta), "=", round(coef[2,1],2))) +
geom_text(x=3000, y=800, label=paste("p-value", "=", round(coef[2,4],6)))+
xlab("January Incidence Rates")+ylab("August Incidence Rates")
## `geom_smooth()` using formula 'y ~ x'
CityHispanic.cont <- as.numeric(relevel(CityHispanic.cat, "Low"))
coef <- summary(lm(Aug.IR~CityHispanic.cont))$coef
ggplot(d.analysis, aes(x=CityHispanic.cont, y=Aug.IR)) +
geom_point(color="blue", size=2) +
geom_smooth(method=lm, color="black") +
geom_text(x=2.5, y=900, label=paste(expression(beta), "=", round(coef[2,1],2))) +
geom_text(x=2.5, y=800, label=paste("p-value", "=", round(coef[2,4],2))) +
scale_x_continuous(name ="Proportion of Hispanic Individuals", breaks=c(1,2,3), labels=c("Low", "Med", "High"))+
ylab("August Incidence Rates")
## `geom_smooth()` using formula 'y ~ x'
coef <- summary(lm(Aug.IR~as.numeric(relevel(HouseholdIncome.cat, "Low"))))$coef
ggplot(d.analysis, aes(x=as.numeric(relevel(HouseholdIncome.cat, "Low")), y=Aug.IR)) +
geom_point(color="blue", size=2) +
geom_smooth(method=lm, color="black") +
geom_text(x=2.5, y=900, label=paste(expression(beta), "=", round(coef[2,1],2))) +
geom_text(x=2.5, y=800, label=paste("p-value", "=", round(coef[2,4],2))) +
scale_x_continuous(name ="Household Income", breaks=c(1,2,3), labels=c("Low", "Med", "High"))+
ylab("August Incidence Rates")
## `geom_smooth()` using formula 'y ~ x'
coef <- summary(lm(Aug.IR~Aug.Vax.percent))$coef
ggplot(d.analysis, aes(x=Aug.Vax.percent, y=Aug.IR)) +
geom_point(color="blue", size=2) +
geom_smooth(method=lm, color="black") +
geom_text(x=80, y=900, label=paste(expression(beta), "=", round(coef[2,1],2))) +
geom_text(x=80, y=800, label=paste("p-value", "=", round(coef[2,4],2)))+
xlab("August Vaccination Percent")+ylab("August Incidence Rates")
## `geom_smooth()` using formula 'y ~ x'
reg.Jan.noVax <- lm(Jan.IR ~ HouseholdIncome.cat+CityRaceEthnicty.m+CityAge.m)
kable(summary(reg.Jan.noVax)$coef, digits=3, format = "simple", align = 'c', caption = "January Regression")
| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 1087.329 | 431.053 | 2.522 | 0.013 |
| HouseholdIncome.catHigh | -983.811 | 171.194 | -5.747 | 0.000 |
| HouseholdIncome.catLow | 783.304 | 163.598 | 4.788 | 0.000 |
| CityRaceEthnicty.mCityHispanic.cat | 97.188 | 81.787 | 1.188 | 0.237 |
| CityRaceEthnicty.mCityBlack.cat | 65.163 | 165.618 | 0.393 | 0.695 |
| CityRaceEthnicty.mCityAsian.cat | 67.657 | 73.121 | 0.925 | 0.357 |
| CityAge.mCityAgeOld.cat | -107.726 | 146.655 | -0.735 | 0.464 |
| CityAge.mCityAgeYoung.cat | 444.903 | 164.853 | 2.699 | 0.008 |
#stargazer(reg.Jan.noVax, type = 'html', ci=TRUE, ci.level=0.95, single.row = T)
\(R^2 =\) 0.62
kable(anova(reg.Jan.noVax), digits=3, format = "simple", align = 'c', caption = "January Regression ANOVA")
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| HouseholdIncome.cat | 2 | 60973275.2 | 30486637.6 | 87.403 | 0.000 |
| CityRaceEthnicty.m | 3 | 542565.2 | 180855.1 | 0.518 | 0.670 |
| CityAge.m | 2 | 3052626.7 | 1526313.3 | 4.376 | 0.015 |
| Residuals | 113 | 39415190.7 | 348807.0 | NA | NA |
reg.Aug.vax <- lm(Aug.IR ~ Aug.Vax.cat)
kable(summary(reg.Aug.vax)$coef, digits=3, format = "simple", align = 'c', caption = "August Regression")
| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 450.118 | 21.753 | 20.692 | 0.000 |
| Aug.Vax.catHigh | -109.075 | 37.886 | -2.879 | 0.005 |
| Aug.Vax.catLow | 54.182 | 37.886 | 1.430 | 0.155 |
\(R^2 =\) 0.11
kable(anova(reg.Aug.vax), digits=3, format = "simple", align = 'c', caption = "August Regression ANOVA")
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| Aug.Vax.cat | 2 | 422581 | 211290.5 | 7.32 | 0.001 |
| Residuals | 118 | 3406082 | 28865.1 | NA | NA |
reg.Aug <- lm(Aug.IR ~ Aug.Vax.cat + Jan.IR.cat + HouseholdIncome.cat+CityRaceEthnicty.m+CityAge.m)
kable(summary(reg.Aug)$coef, digits=3, format = "simple", align = 'c', caption = "August Regression")
| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 298.772 | 114.632 | 2.606 | 0.010 |
| Aug.Vax.catHigh | -53.373 | 42.644 | -1.252 | 0.213 |
| Aug.Vax.catLow | 43.507 | 41.916 | 1.038 | 0.302 |
| Jan.IR.catHigh | -37.110 | 47.008 | -0.789 | 0.432 |
| Jan.IR.catLow | -218.582 | 49.835 | -4.386 | 0.000 |
| HouseholdIncome.catHigh | 59.182 | 48.820 | 1.212 | 0.228 |
| HouseholdIncome.catLow | -67.747 | 48.738 | -1.390 | 0.167 |
| CityRaceEthnicty.mCityHispanic.cat | 73.615 | 25.706 | 2.864 | 0.005 |
| CityRaceEthnicty.mCityBlack.cat | 88.676 | 44.028 | 2.014 | 0.046 |
| CityRaceEthnicty.mCityAsian.cat | -5.201 | 20.844 | -0.250 | 0.803 |
| CityAge.mCityAgeOld.cat | -70.434 | 39.484 | -1.784 | 0.077 |
| CityAge.mCityAgeYoung.cat | 58.694 | 45.717 | 1.284 | 0.202 |
\(R^2 =\) 0.31
kable(anova(reg.Aug), digits=3, format = "simple", align = 'c', caption = "August Regression ANOVA")
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| Aug.Vax.cat | 2 | 422581.0 | 211290.50 | 8.714 | 0.000 |
| Jan.IR.cat | 2 | 233137.6 | 116568.79 | 4.808 | 0.010 |
| HouseholdIncome.cat | 2 | 98384.8 | 49192.40 | 2.029 | 0.136 |
| CityRaceEthnicty.m | 3 | 287587.3 | 95862.42 | 3.954 | 0.010 |
| CityAge.m | 2 | 144166.4 | 72083.21 | 2.973 | 0.055 |
| Residuals | 109 | 2642806.1 | 24245.93 | NA | NA |
reg.Aug.vax <- lm(Aug.Vax.percent ~ Jan.IR.cat+HouseholdIncome.cat+CityRaceEthnicty.m+CityAge.m)
kable(summary(reg.Aug.vax)$coef, digits=3, format = "simple", align = 'c', caption = "August Vaccination Regression")
| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 59.663 | 6.915 | 8.627 | 0.000 |
| Jan.IR.catHigh | -1.763 | 2.791 | -0.632 | 0.529 |
| Jan.IR.catLow | 0.389 | 2.866 | 0.136 | 0.892 |
| HouseholdIncome.catHigh | 3.236 | 2.958 | 1.094 | 0.276 |
| HouseholdIncome.catLow | -0.898 | 2.940 | -0.305 | 0.761 |
| CityRaceEthnicty.mCityHispanic.cat | 3.718 | 1.521 | 2.445 | 0.016 |
| CityRaceEthnicty.mCityBlack.cat | 0.993 | 2.651 | 0.375 | 0.709 |
| CityRaceEthnicty.mCityAsian.cat | 3.690 | 1.176 | 3.137 | 0.002 |
| CityAge.mCityAgeOld.cat | -2.160 | 2.381 | -0.907 | 0.366 |
| CityAge.mCityAgeYoung.cat | -6.132 | 2.681 | -2.287 | 0.024 |
kable(anova(reg.Aug.vax), digits=3, format = "simple", align = 'c', caption = "August Vaccination Regression ANOVA")
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| Jan.IR.cat | 2 | 1976.111 | 988.056 | 11.094 | 0.000 |
| HouseholdIncome.cat | 2 | 1059.438 | 529.719 | 5.948 | 0.004 |
| CityRaceEthnicty.m | 3 | 1696.679 | 565.560 | 6.350 | 0.001 |
| CityAge.m | 2 | 488.550 | 244.275 | 2.743 | 0.069 |
| Residuals | 111 | 9885.456 | 89.058 | NA | NA |
reg.Aug.vax <- lm(Aug.Vax.percent ~ Jan.IR.cat+HouseholdIncome.cat+CityRaceEthnicty.m+CityAge.m)
reg.Aug.vax <- multinom(Aug.Vax.cat ~ Jan.IR.cat+HouseholdIncome.cat+CityRaceEthnicty.m+CityAge.m)
## # weights: 33 (20 variable)
## initial value 132.932087
## iter 10 value 78.668921
## iter 20 value 76.286634
## iter 30 value 76.182929
## iter 40 value 76.181907
## final value 76.181903
## converged
coef.high <- summary(reg.Aug.vax)$coefficients["High",]
OR.high <- exp(coef.high)
se.high <- summary(reg.Aug.vax)$standard.errors["High",]
z.high <- coef.high/se.high
p.value.high <- (1 - pnorm(abs(z.high), 0, 1)) * 2
r.high <- data.frame(OR.high, coef.high, se.high, z.high, p.value.high)
kable(r.high, digits=3, format = "simple", align = 'c', caption = "August Vaccination Regression for High Vax Rates")
| OR.high | coef.high | se.high | z.high | p.value.high | |
|---|---|---|---|---|---|
| (Intercept) | 459743.494 | 13.038 | 1.130 | 11.538 | 0.000 |
| Jan.IR.catHigh | 0.000 | -16.605 | 0.000 | -12741324.270 | 0.000 |
| Jan.IR.catLow | 6.620 | 1.890 | 0.742 | 2.548 | 0.011 |
| HouseholdIncome.catHigh | 0.783 | -0.245 | 0.834 | -0.294 | 0.769 |
| HouseholdIncome.catLow | 2.879 | 1.058 | 1.343 | 0.787 | 0.431 |
| CityRaceEthnicty.mCityHispanic.cat | 1.589 | 0.463 | 0.413 | 1.122 | 0.262 |
| CityRaceEthnicty.mCityBlack.cat | 0.000 | -15.583 | 1.130 | -13.790 | 0.000 |
| CityRaceEthnicty.mCityAsian.cat | 1.028 | 0.028 | 0.392 | 0.071 | 0.944 |
| CityAge.mCityAgeOld.cat | 1.604 | 0.473 | 0.682 | 0.693 | 0.488 |
| CityAge.mCityAgeYoung.cat | 1.068 | 0.065 | 1.265 | 0.052 | 0.959 |
coef.low <- summary(reg.Aug.vax)$coefficients["Low",]
OR.low <- exp(coef.low)
se.low <- summary(reg.Aug.vax)$standard.errors["Low",]
z.low <- coef.low/se.low
p.value.low <- (1 - pnorm(abs(z.low), 0, 1)) * 2
r.low <- data.frame(OR.low, coef.low, se.low, z.low, p.value.low)
kable(r.low, digits=3, format = "simple", align = 'c', caption = "August Vaccination Regression for Low Vax Rates")
| OR.low | coef.low | se.low | z.low | p.value.low | |
|---|---|---|---|---|---|
| (Intercept) | 3.841 | 1.346 | 2.349 | 0.573 | 0.567 |
| Jan.IR.catHigh | 0.791 | -0.234 | 0.942 | -0.249 | 0.803 |
| Jan.IR.catLow | 4.193 | 1.433 | 1.110 | 1.291 | 0.197 |
| HouseholdIncome.catHigh | 0.533 | -0.629 | 1.031 | -0.610 | 0.542 |
| HouseholdIncome.catLow | 1.750 | 0.560 | 0.934 | 0.600 | 0.549 |
| CityRaceEthnicty.mCityHispanic.cat | 0.377 | -0.976 | 0.574 | -1.701 | 0.089 |
| CityRaceEthnicty.mCityBlack.cat | 0.861 | -0.149 | 0.708 | -0.211 | 0.833 |
| CityRaceEthnicty.mCityAsian.cat | 0.081 | -2.516 | 0.799 | -3.151 | 0.002 |
| CityAge.mCityAgeOld.cat | 2.407 | 0.878 | 0.853 | 1.029 | 0.303 |
| CityAge.mCityAgeYoung.cat | 3.720 | 1.314 | 0.807 | 1.629 | 0.103 |
Full model fit by date
{r} adj_case_14day_rate.x ~ Vax.percent + Jan.IR.cat + HouseholdIncome.cat+CityHispanic.cat+CityBlack.cat+ CityAsian.cat+ CityAgeOld.cat+ CityAgeYoung.cat
No vaccination lag time in this example
Modeled as a continuous variable. Plots include dates > 04/01/2021
d <- readRDS("data/analysis_data.rds")
lac_summarytable_dash <- read.csv("data/LA_County_Covid19_CSA_14day_case_death_table_Clean.csv")
# data clean ----
dout <- d %>%
mutate(Date.x = ymd(Date.x),
Vax.percent = dose1_all_c_prcent,
Jan.IR = as.numeric(adj_case_14day_rate.y),
Jan.IR.cat = relevel(as.factor(ifelse(Jan.IR > quantile(Jan.IR, probs=.75), "High",
ifelse(Jan.IR < quantile(Jan.IR, probs=.25), "Low",
"Med"))), ref="Low")) %>%
dplyr::select(Date.x, city, adj_case_14day_rate.x,
Vax.percent, Jan.IR.cat,
HouseholdIncome.cat, CityHispanic.cat, CityBlack.cat,
CityAsian.cat, CityAgeOld.cat, CityAgeYoung.cat)
test <- dout %>%
# dplyr::filter(Date.x > "2021-03-01") %>%
tidyr::nest(data = -Date.x) %>%
dplyr::mutate(fit = purrr::map(data, ~lm(adj_case_14day_rate.x ~ Vax.percent + Jan.IR.cat + HouseholdIncome.cat+CityHispanic.cat+CityBlack.cat+ CityAsian.cat+ CityAgeOld.cat+ CityAgeYoung.cat, data = .x)),
tidied = purrr::map(fit, ~broom::tidy(.x), quality = purrr::map(fit, ~glance(.x)))) %>%
dplyr::select(-data, -fit) %>%
tidyr::unnest(tidied) %>%
# dplyr::filter(grepl("Vax", term)) %>%
arrange(Date.x)
# Vaccination rates ----
# after april 2021 (vaccination rates too low prior)
p1 <- ggplot(test %>% filter(grepl("Vax.percent", term), Date.x > "2021-04-01"), aes(Date.x, estimate)) +
geom_point() +
geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error)) +
geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.8))
# geom_line(data = lac_summarytable_dash %>% filter(city == "Los Angeles County"), aes(Date, adj_case_14day_rate))
p2 <- ggplot(lac_summarytable_dash %>%
filter(city == "Los Angeles County",
Date > "2021-04-01") %>%
mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
geom_line()
plot_grid(p1, p2, ncol = 1, align = 'v')
p1 <- ggplot(test %>% filter(term %in% c("CityAsian.catHigh", "CityBlack.catHigh", "CityHispanic.catHigh")), aes(Date.x, estimate)) +
geom_point(aes(colour = term)) +
geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error, colour = term)) +
geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.8))
p2 <- ggplot(lac_summarytable_dash %>%
filter(city == "Los Angeles County") %>%
mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
geom_line()
plot_grid(p1, p2, ncol = 1, align = 'v')
p1 <- ggplot(test %>% filter(term %in% c("HouseholdIncome.catHigh", "HouseholdIncome.catLow")), aes(Date.x, estimate)) +
geom_point(aes(colour = term)) +
geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error, colour = term)) +
geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.2))
p2 <- ggplot(lac_summarytable_dash %>%
filter(city == "Los Angeles County") %>%
mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
geom_line()
plot_grid(p1, p2, ncol = 1, align = 'v')
## Warning: Removed 1 row(s) containing missing values (geom_path).
p1 <- ggplot(test %>% filter(term %in% c("CityAgeOld.catHigh", "CityAgeYoung.catHigh")), aes(Date.x, estimate)) +
geom_point(aes(colour = term)) +
geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error, colour = term)) +
geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.8))
p2 <- ggplot(lac_summarytable_dash %>%
filter(city == "Los Angeles County") %>%
mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
geom_line()
plot_grid(p1, p2, ncol = 1, align = 'v')
## Warning: Removed 1 row(s) containing missing values (geom_path).
p1 <- ggplot(test %>% filter(term %in% c("Jan.IR.catHigh", "Jan.IR.catMed")), aes(Date.x, estimate)) +
geom_point(aes(colour = term)) +
geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error, colour = term)) +
geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.8))
# geom_line(data = lac_summarytable_dash %>% filter(city == "Los Angeles County"), aes(Date, adj_case_14day_rate))
# ggsave(p1, width = 10, height = 8, file = "~/Desktop/test1.png")
p2 <- ggplot(lac_summarytable_dash %>%
filter(city == "Los Angeles County") %>%
mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
geom_line()
plot_grid(p1, p2, ncol = 1, align = 'v')
## Warning: Removed 1 row(s) containing missing values (geom_path).
Full model fit by date
{r} adj_case_14day_rate.x ~ Vax.percent + Jan.IR.cat + HouseholdIncome.cat+CityHispanic.cat+CityBlack.cat+ CityAsian.cat+ CityAgeOld.cat+ CityAgeYoung.cat
Introduce 14 day lag between vaccination rates and COVID incidence rates
Modeled as a continuous variable. Plots include dates > 04/01/2021
d <- readRDS("data/analysis_data_14vax_lag.rds")
lac_summarytable_dash <- read.csv("data/LA_County_Covid19_CSA_14day_case_death_table_Clean.csv")
# data clean ----
dout <- d %>%
mutate(Date.x = ymd(Date.x),
Vax.percent = dose1_all_c_prcent,
Jan.IR = as.numeric(adj_case_14day_rate.y),
Jan.IR.cat = relevel(as.factor(ifelse(Jan.IR > quantile(Jan.IR, probs=.75), "High",
ifelse(Jan.IR < quantile(Jan.IR, probs=.25), "Low",
"Med"))), ref="Low")) %>%
dplyr::select(Date.x, city, adj_case_14day_rate.x,
Vax.percent, Jan.IR.cat,
HouseholdIncome.cat, CityHispanic.cat, CityBlack.cat,
CityAsian.cat, CityAgeOld.cat, CityAgeYoung.cat)
test <- dout %>%
# dplyr::filter(Date.x > "2021-03-01") %>%
tidyr::nest(data = -Date.x) %>%
dplyr::mutate(fit = purrr::map(data, ~lm(adj_case_14day_rate.x ~ Vax.percent + Jan.IR.cat + HouseholdIncome.cat+CityHispanic.cat+CityBlack.cat+ CityAsian.cat+ CityAgeOld.cat+ CityAgeYoung.cat, data = .x)),
tidied = purrr::map(fit, ~broom::tidy(.x), quality = purrr::map(fit, ~glance(.x)))) %>%
dplyr::select(-data, -fit) %>%
tidyr::unnest(tidied) %>%
# dplyr::filter(grepl("Vax", term)) %>%
arrange(Date.x)
# Vaccination rates ----
# after april 2021 (vaccination rates too low prior)
p1 <- ggplot(test %>% filter(grepl("Vax.percent", term), Date.x > "2021-04-01"), aes(Date.x, estimate)) +
geom_point() +
geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error)) +
geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.8))
# geom_line(data = lac_summarytable_dash %>% filter(city == "Los Angeles County"), aes(Date, adj_case_14day_rate))
p2 <- ggplot(lac_summarytable_dash %>%
filter(city == "Los Angeles County",
Date > "2021-04-01") %>%
mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
geom_line()
plot_grid(p1, p2, ncol = 1, align = 'v')
p1 <- ggplot(test %>% filter(term %in% c("CityAsian.catHigh", "CityBlack.catHigh", "CityHispanic.catHigh")), aes(Date.x, estimate)) +
geom_point(aes(colour = term)) +
geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error, colour = term)) +
geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.8))
p2 <- ggplot(lac_summarytable_dash %>%
filter(city == "Los Angeles County") %>%
mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
geom_line()
plot_grid(p1, p2, ncol = 1, align = 'v')
p1 <- ggplot(test %>% filter(term %in% c("HouseholdIncome.catHigh", "HouseholdIncome.catLow")), aes(Date.x, estimate)) +
geom_point(aes(colour = term)) +
geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error, colour = term)) +
geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.2))
p2 <- ggplot(lac_summarytable_dash %>%
filter(city == "Los Angeles County") %>%
mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
geom_line()
plot_grid(p1, p2, ncol = 1, align = 'v')
## Warning: Removed 1 row(s) containing missing values (geom_path).
p1 <- ggplot(test %>% filter(term %in% c("CityAgeOld.catHigh", "CityAgeYoung.catHigh")), aes(Date.x, estimate)) +
geom_point(aes(colour = term)) +
geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error, colour = term)) +
geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.8))
p2 <- ggplot(lac_summarytable_dash %>%
filter(city == "Los Angeles County") %>%
mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
geom_line()
plot_grid(p1, p2, ncol = 1, align = 'v')
## Warning: Removed 1 row(s) containing missing values (geom_path).
p1 <- ggplot(test %>% filter(term %in% c("Jan.IR.catHigh", "Jan.IR.catMed")), aes(Date.x, estimate)) +
geom_point(aes(colour = term)) +
geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error, colour = term)) +
geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.8))
# geom_line(data = lac_summarytable_dash %>% filter(city == "Los Angeles County"), aes(Date, adj_case_14day_rate))
# ggsave(p1, width = 10, height = 8, file = "~/Desktop/test1.png")
p2 <- ggplot(lac_summarytable_dash %>%
filter(city == "Los Angeles County") %>%
mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
geom_line()
plot_grid(p1, p2, ncol = 1, align = 'v')
## Warning: Removed 1 row(s) containing missing values (geom_path).